gbsv MC2¶
Day 1 - Auto-Correlation¶
Overarching Task: Choose a country and theme all the MC tasks to this county. The use cases, problem statements, signals and images should be related to your country of choice. It can be the same country as in MC1, but you may also choose a different one.
Day 1 Task Passend zu deinem Land: Definiere einen Anwendungsfall, um zu erkennen, ob ein 1D Signal wiederkehrende Muster enthält. Suche ein passendes 1D Signal (Audio, Zeitreihe, Vitalparameter, ...), welches wiederkehrende Muster enthält, um in den Folgetagen Auto-Korrelation anzuwenden.
import time
start_time = time.time()
import cv2 as cv
import librosa
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import sounddevice as sd
from ipywidgets import interact, widgets
from matplotlib.widgets import Slider, Button, TextBox
from matplotlib.patches import Rectangle
from scipy.io.wavfile import write
from scipy.signal import correlate
from IPython.display import Markdown, Audio
from rich import print
from rich import pretty
pretty.install()
audio, sample_rate = librosa.load("sound/bell.wav")
print(f"{audio=}")
print(f"{audio.dtype=}")
print(f"{audio.shape=}")
print(f"{np.max(audio)=}")
print(f"{np.min(audio)=}")
display(Markdown(f"**ring the Bell**"))
display(Audio(audio, rate=sample_rate))
audio=array([2.6749703e-35, 8.3550063e-35, 1.1725352e-34, ..., 6.8832678e-04, 7.8442483e-04, 0.0000000e+00], dtype=float32)
audio.dtype=dtype('float32')
audio.shape=(1641260,)
np.max(audio)=np.float32(0.86444175)
np.min(audio)=np.float32(-0.8645316)
ring the Bell
The chosen use case is to estimate the weight/size of a Swiss cowbell based on its frequency characteristics using auto-correlation. Cowbells produce distinct resonant frequencies determined by their size, material, and weight. Identifying recurring patterns in their acoustic signal can help derive these physical properties. i chose this sound because cowbells are an iconic symbol of Swiss craftsmanship and tradition, and we can always hear them on our hiking trips on weekends. The objective of this experiment is to demonstrate whether auto-correlation can reliably detect these recurring patterns and link them to the bell’s weight. The chosen signal is a recording of a Swiss cowbell.
Day 2 Task¶
Analysiere mittels Auto-Korrelation die wiederkehrenden Muster innerhalb deines 1D Signals.
times = np.arange(0, len(audio)) / sample_rate
plt.figure(figsize=(12, 4))
plt.plot(times, audio)
plt.title("Original Bell Sound")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.show()
First i will analyze the signal in the time domain, the autocorrelation, and the frequency spectrum for two cases: one signal includes the initial hit of the bell, while the other excludes it. This helps us determine how the initial strike influences the periodicity and dominant frequencies.
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from scipy.signal import correlate, find_peaks
from scipy.fftpack import fft
from IPython.display import Markdown, Audio
def trim_signal(signal, sample_rate, start_time, end_time):
start_sample = int(start_time * sample_rate)
end_sample = int(end_time * sample_rate)
signal = signal[start_sample:end_sample]
return signal
signal = trim_signal(audio, sample_rate, 1.35, 1.5)
write('sound/bell_no_hit.wav', sample_rate, signal)
display(Markdown(f"## **this is only the ringing, without the initial hit**"))
display(Audio(signal, rate=sample_rate))
times = np.arange(0, len(signal)) / sample_rate
plt.figure(figsize=(12, 4))
plt.plot(times, signal)
plt.title("this is only the ringing, without the initial hit")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.show()
signal_hit = trim_signal(audio, sample_rate, 1.05, 1.5)
write('sound/bell_with_hit.wav', sample_rate, signal)
display(Markdown(f"## **this is the ringing, with the initial hit**"))
display(Audio(signal_hit, rate=sample_rate))
times = np.arange(0, len(signal_hit)) / sample_rate
plt.figure(figsize=(12, 4))
plt.plot(times, signal_hit)
plt.title("this is the ringing, with the initial hit")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.show()
def plot_signal_analysis(signal1, signal2=None, sample_rate=None, labels=("Signal 1", "Signal 2")):
def compute_autocorrelation(signal):
autocorr = correlate(signal, signal, mode='full')
lags = np.arange(-len(signal) + 1, len(signal))
return lags, autocorr
def calculate_kpis(signal):
mean = np.mean(signal)
std = np.std(signal)
var = np.var(signal)
return {"Mean": mean, "Std Dev": std, "Variance": var}
def compute_fft(signal, sample_rate):
N = len(signal)
freqs = np.fft.fftfreq(N, d=1/sample_rate)
fft_values = np.abs(fft(signal))[:N//2]
return freqs[:N//2], fft_values
signals = [signal1] if signal2 is None else [signal1, signal2]
num_signals = len(signals)
fig, axes = plt.subplots(3, num_signals, figsize=(7 * num_signals, 12))
axes = axes if num_signals > 1 else [axes]
for idx, signal in enumerate(signals):
# Time Domain
times = np.arange(len(signal)) / sample_rate
axes[0][idx].plot(times, signal)
axes[0][idx].set_title(f"{labels[idx]} - Time Domain")
axes[0][idx].set_xlabel("Time (s)")
axes[0][idx].set_ylabel("Amplitude")
axes[0][idx].grid()
# Auto-Correlation
lags, autocorr = compute_autocorrelation(signal)
peaks, _ = find_peaks(autocorr[len(autocorr)//2:], distance=50)
peak_lags = lags[len(autocorr)//2:][peaks]
axes[1][idx].plot(lags, autocorr)
axes[1][idx].scatter(peak_lags, autocorr[len(autocorr)//2:][peaks], color='red', label='Peaks')
axes[1][idx].set_title(f"{labels[idx]} - Auto-Correlation")
axes[1][idx].set_xlabel("Lag")
axes[1][idx].set_ylabel("Correlation")
axes[1][idx].legend()
axes[1][idx].grid()
# Frequency Spectrum (FFT)
freqs, fft_values = compute_fft(signal, sample_rate)
axes[2][idx].plot(freqs, fft_values)
axes[2][idx].set_title(f"{labels[idx]} - Frequency Spectrum")
axes[2][idx].set_xlabel("Frequency (Hz)")
axes[2][idx].set_ylabel("Amplitude")
axes[2][idx].grid()
# Identify peaks in the FFT values
peaks, properties = find_peaks(fft_values, height=100)
dominant_frequencies = freqs[peaks]
# Display the frequencies
print(f"{labels[idx]} Dominant Frequencies (Hz):", dominant_frequencies)
# KPIs
kpis = calculate_kpis(signal)
print(f"KPIs for {labels[idx]}: {kpis}")
plt.tight_layout()
plt.show()
plot_signal_analysis(signal1=signal, signal2=signal_hit, sample_rate=sample_rate, labels=("Ringing without hit", "Ringing with hit"))
Ringing without hit Dominant Frequencies (Hz): [1379.79141475]
KPIs for Ringing without hit: {'Mean': np.float32(-0.0001495729), 'Std Dev': np.float32(0.06533235), 'Variance': np.float32(0.004268316)}
Ringing with hit Dominant Frequencies (Hz): [ 328.87231684 1379.93046458 2095.44996473]
KPIs for Ringing with hit: {'Mean': np.float32(-6.990738e-05), 'Std Dev': np.float32(0.15481316), 'Variance': np.float32(0.023967113)}
The signal with the initial hit demonstrates a sharper amplitude decay and additional low-frequency content, which is evident in the frequency spectrum where lower frequencies dominate the early part of the signal. The autocorrelation plot for the signal without the hit displays a smoother and more consistent decay, highlighting stronger periodicity and indicating that the ringing alone is more stable. This conclusion is supported by the calculated KPIs: the signal without the hit has a lower variance (0.0043 compared to 0.0239) and standard deviation (0.065 compared to 0.154), which indicate more predictable and stable oscillations. In contrast, the initial strike introduces transient energy and broadens the frequency spectrum, adding harmonics but reducing the overall stability of the periodic pattern, as reflected in the less uniform autocorrelation decay.
with the next plots i will focuse on zoomed-in regions of the ringing signal, comparing the raw signal to aligned sine waves of the dominant frequencies extracted from the FFT. This demonstrates how well these frequencies represent the periodic structure of the signal.
from scipy.signal import correlate
zoom_stage1_start, zoom_stage1_end = 0.02, 0.04
zoom_stage2_start, zoom_stage2_end = 0.09, 0.11
times = np.arange(0, len(signal_hit)) / sample_rate
zoom_times_stage1 = times[int(zoom_stage1_start * sample_rate):int(zoom_stage1_end * sample_rate)]
zoom_signal_stage1 = signal_hit[int(zoom_stage1_start * sample_rate):int(zoom_stage1_end * sample_rate)]
zoom_times_stage2 = times[int(zoom_stage2_start * sample_rate):int(zoom_stage2_end * sample_rate)]
zoom_signal_stage2 = signal_hit[int(zoom_stage2_start * sample_rate):int(zoom_stage2_end * sample_rate)]
dominant_frequencies = [328.87231684, 1379.93046458, 2095.44996473]
def calculate_phase_offset(signal, sine_wave):
correlation = correlate(signal, sine_wave, mode='full')
lag = np.argmax(correlation) - len(sine_wave) + 1
return lag
def generate_aligned_sine_wave(frequency, time, signal, amplitude=0.2):
sine_wave = amplitude * np.sin(2 * np.pi * frequency * time)
lag = calculate_phase_offset(signal, sine_wave)
aligned_sine_wave = np.roll(sine_wave, lag)
return aligned_sine_wave
# Stage 2: Zoom Region 1
plt.figure(figsize=(12, 12))
plt.subplot(3, 1, 1)
plt.plot(times, signal_hit, label="Signal")
plt.axvspan(zoom_stage1_start, zoom_stage1_end, color='red', alpha=0.3, label="Zoom Region 1")
plt.axvspan(zoom_stage2_start, zoom_stage2_end, color='green', alpha=0.3, label="Zoom Region 2")
plt.title("Full Signal with Zoom Regions Highlighted")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.ylim(-1, 1)
plt.legend()
plt.grid()
plt.subplot(3, 1, 2)
aligned_sine_wave = generate_aligned_sine_wave(
dominant_frequencies[0], zoom_times_stage1[:len(zoom_signal_stage1)], zoom_signal_stage1, 0.75
)
plt.plot(zoom_times_stage1, aligned_sine_wave, color='red', label=f"Sine {dominant_frequencies[0]:.1f} Hz (Aligned)")
aligned_sine_wave = generate_aligned_sine_wave(
dominant_frequencies[1], zoom_times_stage1[:len(zoom_signal_stage1)], zoom_signal_stage1
)
plt.plot(zoom_times_stage1, aligned_sine_wave, color='orange', label=f"Sine {dominant_frequencies[1]:.1f} Hz (Aligned)")
plt.plot(zoom_times_stage1, zoom_signal_stage1, color='blue', label="Signal")
plt.title("Zoom Region 1 with Aligned Dominant Frequencies")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.ylim(-1, 1)
plt.legend()
plt.grid()
# Stage 3: Zoom Region 2
plt.subplot(3, 1, 3)
aligned_sine_wave = generate_aligned_sine_wave(
dominant_frequencies[1], zoom_times_stage2[:len(zoom_signal_stage2)], zoom_signal_stage2
)
plt.plot(zoom_times_stage2, aligned_sine_wave, color='orange', label=f"Sine {dominant_frequencies[1]:.1f} Hz (Aligned)")
aligned_sine_wave = generate_aligned_sine_wave(
dominant_frequencies[2], zoom_times_stage2[:len(zoom_signal_stage2)], zoom_signal_stage2
)
plt.plot(zoom_times_stage2, aligned_sine_wave-0.5, color='green', label=f"Sine {dominant_frequencies[2]:.1f} Hz (Aligned)")
plt.plot(zoom_times_stage2, zoom_signal_stage2, color='blue', label="Signal")
plt.title("Zoom Region 2 with Aligned Dominant Frequencies")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.ylim(-1, 1)
plt.legend()
plt.grid()
plt.tight_layout()
plt.show()
Observation: In the first zoomed region, the sine wave at 328.9 Hz closely follows the signal's overall oscillations, showing its dominance in the lower-frequency range. This can only be observed right after the initial hit and fades the fastest. In the second zoomed region, the sine wave at 1380 Hz captures finer details, confirming its contribution to the higher-frequency components. These aligned sine waves visually highlight the signal's harmonic structure. Note that 1380 and 2095 Hz Frequencies are alredy present in the beginning but fade much slower. The 2095 Hz frequencie was plotted below the other ones to prevent overcrouding of the plot, but is is still interesting to observe its contribution.